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Abstract 

We present a significant improvement to a time-dependent WKB (TDWKB) formulation developed 
by Boiron and Lombardi [JCP 108, 3431 (1998)] in which the TDWKB equations are solved along 
classical trajectories that propagate in the complex plane. Boiron and Lombardi showed that 
the method gives very good agreement with the exact quantum mechanical result as long as the 
wavefunction does not exhibit interference effects such as oscillations and nodes. In this paper we 
show that this limitation can be overcome by superposing the contributions of crossing trajectories. 
We also demonstrate that the approximation improves when incorporating higher order terms in 
the expansion. These improvements could make the TDWKB formulation a competitive alternative 
to current time-dependent semiclassical methods. 
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I. INTRODUCTION 



The difficulty in performing quantum mechanical calculations of multi-dimensional systems 
has stimulated an intensive and ongoing effort in the last three decades to develop numerical 
tools based on semiclassical mechanics. In this context, we refer to semiclassical mechanics 
as the derivation of a quantum mechanical wavefunction or propagator via propagation of 
classical (or classical-like) trajectories. From a physical point of view, semiclassical methods 
try to evade the non-locality imbedded in quantum mechanics. Mathematically speaking, 
semiclassical methods aim at casting the time- dependent Schrodinger equation (TDSE), 
which is a PDE, in terms of ODEs related to classical equations of motion. This trans- 
formation has significant computational advantages that can ease the inherent difficulty of 
multi-dimensional quantum calculations. 

The WKB method[l|, O, can be considered as the ffist of the semiclassical methods. Its 
date of birth almost coincides with the publication of the Schrodinger equation in 1926, and 
virtually every standard text book in quantum mechanics has a description of the method. 
The basic idea of the WKB method is to recast the wavefunction as the exponential of a 
function and then replace the exponent with a power series in h. The WKB method is 
ordinarily applied to the time-independent Schrodinger equation and provides for a good 
approximation to the eigenstates as long as one is not too near a classical turning point. It 
is only natural that as part of the effort to develop time-dependent semiclassical methods. 



a time- dependent version of the WKB method wou 



has been done in this direction Q, B, y, 0, y, y, q, q, q, q, y, h, q. A decade ago. 



3e explored. Surprisingly little work 



Boiron and Lombardi 



171] developed a complex trajectory version of time- dependent WKB 



(TDWKB), which we refer to as CTDWKB. In conventional WKB the leading order term 
in the phase of the wave function is taken to be 0{h~^) and the leading order term in the 
amplitude is taken to be 0{hP). In contrast, the CTDWKB formulation treats the amplitude 
and phase on an equal footing. The price to pay for this procedure is that the resulting 
classical trajectories propagate in the complex plane. The benefits are that the results are 
superior to standard TDWKB and no singularities are encountered during the integration 
of the equation of motion. 

The CTDWKB equations of motion can be solved analytically and yield the exact wave- 
function for an initial Gaussian wavepacket in a potential with up to quadratic terms. The 
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first-order method was tested numerically by Boiron and Lombardi for scattering of a Gaus- 
sian wavepacket from a potential barrier. They showed that the method produced very 
good results as long as the wavefunction did not exhibit interference effects in the form of 
oscillations or nodes[17]. In this paper we present a simple modification to CTDWKB that 
provides an accurate description of oscillations in the wavefunction. We show that complex 
classical trajectories, similar to real classical trajectories, can cross in configuration space. 
By superposing the contributions from two or more crossing trajectories, interference effects 
are obtained. We take CTDWKB a step further in another direction by showing that the ap- 
proximation generally improves when incorporating additional terms in the series expansion. 
Since the WKB expansion is an asymptotic series, this observation is non-trivial. 

Two other semiclassical formulations that incorporate complex trajectories should be 
mentioned in relation with CTDWKB. The first is the Generalized Gaussian Wavepacket 
Dynamics (GGWPD) developed by Huber, Heller and Littlejohn jisl . One may show 
that for an initial Gaussian wavepacket the equations of motion of GGWPD are de facto 
identical to the equations of the first-order approximation of CTDWKB. However the GG- 
WPD has no generalization to arbitrary initial wavefunctions and no systematic way to 
increase the accuracy of the approximation. On the other hand, in reference p^], Huber and 
Heller appreciate the importance of multiple complex trajectories in obtaining interference 
phenomena. Here we incorporate the idea of crossing complex trajectories into the more 
general CTDWKB formulation. 



The second formulation t 
Complex Action (BOMCA) 



lat is closely related to CTDWKB is Bohmian Mechanics with 
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2l|. CTDWKB and BOMCA begin with the same ansatz 



of substituting of an exponential function ex]i[iS/h] into the TDSE. Similar to CTDWKB, 
the BOMCA formulation incorporates equations of motion that propagate along complex 
trajectories. The first-order equations of motion of BOMCA are identical to the equations of 
first-order CTDWKB. The differences between the two formulations are: (1) The equations 
of motion in BOMCA are for the coefficients of spatial derivatives of the phase. In CTDWKB 
the equations of motion are for the coefficients of an h Taylor expansion of the phase and their 
spatial derivatives. (2) Incorporating higher order terms of the CTDWKB approximation 
does not effect the the results for lower order terms since each equation of motion depends 
only on lower terms of the expansion. This is not the case with BOMCA where each equation 
of motion depends on both lower and higher terms resulting in a backward feedback. (3) A 
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result of the last difference is that in CTDWKB the equations of motion of the trajectories 
remain classical whereas in BOMCA, the inclusion of higher orders of the approximation 
affect the complex trajectories by adding a "quantum force" that yields quantum trajectories. 

This paper is organized as follows. In section [ITl we formulate TDWKB and the CTD- 
WKB. Our derivation is more compact than the Boiron-Lombardi derivation and demon- 
strates how to obtain the equations of motion for higher orders of the expansion in a simple 
manner. In Section UTTl we apply the formulation to a Gaussian initial wavepacket propagat- 
ing in a quartic double- well potential. We demonstrate that superimposing the contributions 
of crossing trajectories leads to interference effects and that incorporating higher order terms 
in the expansion improves the approximation. Section [IV] is a summary and concluding re- 
marks. Following Boiron and Lombardi we will refer to the CTDWKB method in the body 
of the paper as the complex trajectory method (CTM) for short. 

II. FORMULATION 

A. Time-independent vs. Time-dependent WKB 

For simplicity we present the one-dimensional version of the CTM derivation. The general- 
ization to mult i- dimensions can be performed in a straightforward manner. The conventional 
WKB derivation begins by inserting the ansatz 

(2.1) 

into the time-independent Schrodinger equation, where h is Planck's constant divided by 
2n. The end result is 

1 fdSV , ih d'^S ^ 

— — + V(x) — = E, 2.2 

2m\dx J ^ ' 2m dx^ ' ^ ^ 

where m is the mass of the particle, V{x) is the potential energy and E is the eigenvalue. If 
we assume that S{x) can be expanded asymptotically as a polynomial in h 

oo 

S{x) = So{x) + hSi{x) + ffS2{x) + ... = J2 h^Sj{x), (2.3) 

j=0 

then, by substituting the last equation into eq. (l2.2p and equating powers of h, a series of 
coupled ODEs are obtained for the S/s. 



tpi^x) 



exp 
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The time-dependent WKB begins by inserting the ansatz 



m 



ipi^x, t) = exp 



into the time- dependent Schrodinger equation, 



h 



S{x,t) 



ihdti) = -^^xxi^ + v{x, t)ijj, 



The result is the quantum Hamilton- Jacobi equation 



dtS + T^id^S? + V = -^5,,^, 



(2.4) 



(2.5) 



(2.6) 



where the LHS of the equation is in the form of the classical Hamilton- Jacobi equation. 
Equation (12.61) is formally exact since no approximation has been introduced. In TDWKB 
formulation we insert into eq. ( 12. 6p a time-dependent version of eq. ( 12. 3p 



S{x, t) = h? Sj{x, t). 
i=o 



(2.7) 



2m 



(2.8) 



The result is 

oo ^ oo 

j=0 jiJ2=0 ■■" j=0 

By equating terms having the same powers of h we obtain the classical Hamilton- Jacobi 
equation for So{x,t) 

dtSo + ^{dxSoy + V = 0, (2.9) 
and equations of motion for Sn{x, t),n > 1 

n-l 

(2.10) 



dxSo 
m 



dtSn H dxSn — dxxSn—1 ~ ^ dxSj ■ dxSn- 

2m ^-^ 



Conveniently, each equation depends only on lower order terms. The next step in TDWKB 
is to convert eqs. (12.91) and fl2.10p into a set of ODEs by looking at the evolution of 5*0, S\, . . . 
along classical trajectories, as described in the next section. 



B. Integrating along classical trajectories 



As we mentioned earlier, the first term in the U power expansion, 5*0, obeys the Hamilton- 
Jacobi equation (eq. (l2.9l) ). This equation is an alternative formulation of Newton's second 
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law of motion in terms of an action field. The emergence of classical trajectories in the 
TDWKB equations provides the incentive to solve eqs. fl2.9p and (12.101) by integrating along 
such trajectories. 

The link between the Hamilton- Jacobi equation and classical trajectories is demonstrated 
by defining the velocity field 

I +\ — dxSo{x,t) /o 11^ 

v(x,t) = (2.11) 

m 

and considering the trajectories defined by 

^=^(x,t). (2.12) 

By taking the spatial partial derivative of eq. (12.91) . using the definition of the Lagrangian 
time derivative ^ = dt + ^d^, and applying eq. (l2.11l) we obtain the equation of motion for 
the velocity along a trajectory as Newton's second law 

I = (2.13) 

dt m 

Hence, the trajectories defined are simply classical trajectories. 
Inserting eq. fl2.9p in the Lagrangian time derivative of 5*0 yields 

^ = dtSo + ^d.So = \mv' - (2.14) 
dt m 2 

where we recognize the equation of motion for the action along a classical trajectory. Noting 
that f is a mere dummy variable, we summarize the equations of motion for the zeroth order 
term of TDWKB, Sq 

dx dxSo 



dt m 
d{d^So) 



(2.15) 



dt 



d.V, (2.16) 



f^±(9.S.f-V. (2.17) 

We turn to the higher order terms in the series Sn, n > 1. Recognizing the LHS of eqs. fl2.10p 
as the Lagrangian time derivative of Sn, we can write 

n-l 



dSn 

dt 2m '^'^ " ^ 2m 



t; dxxSn~l ~ / dxSj ■ dxSn-j- (2.18) 



These equations do not constitute a closed set of ODEs since they depend on partial deriva- 
tives such as dxxSn- We close the set of equations by deriving equations of motion for the 



partial derivatives on the RHS of eq. fl2.18p {dxxSn-i, dxSj and dxSn-j)- We demonstrate the 
process by deriving equations of motion for 5*1 and 5*2. Inserting n = 1 in eq. fl2.18p yields 

^ = ^dxxSo. (2.19) 
dt 2m ^ ^ ' 

An equation of motion for dxxSo is obtained by taking a second spatial partial derivative of 
eq.dMD, 

dxxtSo + - [dxSo ■ dxxxSo + [dxxSof] + dxxV = 0, (2.20) 
m 



and rewriting it as 



'^-^^^ = --{d..Sof - dxxV. (2.21) 
dt . . m 



This equation is derived in reference 17 by a cumbersome finite difference scheme. It 
is equivalent to eq.(2.9d) of reference ^9] where the equation appears in the context of 
GGWPD. Note that an equation of motion for any order of spatial derivatives of 5*0 can be 
derived in a similar fashion by taking consecutive spatial derivatives of eq. (12.201) and then 
grouping the Lagrangian time derivative terms. Equations fl2.19p and fl2.2ip provide a closed 
set of equations of motion for Si. 
Inserting n = 2 into eq. fl2.18p yields 

^ = :^dxxSi - ^{dxSif. (2.22) 
dt 2m 2m 

The equations of motion for dxSi and dxxSi are obtained by first inserting n = 1 in eq. (l2.10p . 
We then derive two equations by taking a first and a second spatial partial derivative of the 
result. By grouping the Lagrangian time derivatives of dxSi and dxxSi in each of the two 
equations separately we obtain 

— — t: — ^ = '^^dxxxSo dxSi ■ dxxSo, (2.23) 

dt 2m m 

r. = ~ dxxxxSo dxSi ■ dxxxSo dxxSi ■ dxxSo- 

dt 2m m m 

The last equations depend in turn on dxxxSo and dxxxxSo- As mentioned earlier, the equation 
of motion for these terms can be obtained by additional spatial derivatives of eq. fl2.20p . a 
process that yields 

d{dxxxSo) 3 



dt m 



dxxSo ■ dxxxSo — dxxxV, (2.24) 



— [4c}a;3;>S'o ■ dxxxxSo '^(dxxxSo) ] dxxxxV. 

dt m ^ 



Equations (12.210 and fl2.22p - fl2.24p provide a closed set of equations of motion for 5*2. The 
scheme we described for 5*1 and 5*2 can be extended to any of the higher order terms in the 
expansion. Note that incorporating higher order terms Sn in the TDWKB approximation 
does not affect the classical trajectories associated with 5*0, defined by eqs.f l2.l5]) and (12.161) . 
We now turn to the source of the distinction between conventional TDWKB and CTM. 



C. Initial conditions and complex classical trajectories 

In conventional TDWKB the initial wavefunction is "divided" between So{x, 0) and Si{x, 0) 



ip{x, 0) = A{x) exp[i0(x)] = exp 



^So{x,0) + Si{x,0) 



(2.25) 



where A{x) and (f){x) are the initial amplitude and phase respectively, both taken to be real. 
The phase is related to the zero-order term 5*0 and the amplitude to the first-order correction 
term Si according to 

^o(a;,0) = hcpix), Si{x,0) = -i\n[A{x)], (2.26) 

and Sn{x,0) = for > 2. Note that the initial conditions specified by eqs. (12.261) yield 
classical trajectories that propagate on the real axis since 5*0 and its spatial derivatives are 
real quantities (see eqs. fl2.15p and (12.160 ). In contrast, in CTM the amplitude and phase 
are treated on an equal footing with far-reaching consequences. The initial wavefunction is 
specified by So{x,0) 

So{x, 0) = -iMn[i;{x, 0)], S„(x, 0) = 0, n > 1. (2.27) 

Since 5*0 is generally complex and since the initial velocity v{x,0) = dxSo{x,0)/m, the 
trajectories propagate in the complex plane even if the initial positions are on the real axis 
(53[x(0)] = 0). This observation requires us to look at the analytic continuation of the 
wavefunction in the complex plane and find ways to extract the wavefunction on the real 
axis. 



D. Complex root search and superposition 

One of the benefits of conventional TDWKB and CTM compared with BOMCA, is that the 
trajectories obey the classical equations of motion and are independent of the order of the 
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phase expansion we incorporate in the final wavefunction. But the fact still remains that 
for an arbitrary initial position a;(0) G C and an arbitrary final propagation time tf the final 
position x(tf) is complex and yields an "analytically continued" wavefunction at x(tf) 

N 

h 



^lj[x{tf),tf]^exp\lJ2^'^A^(tf),tf]\, (2.28) 



j=0 



where the non-negative integer is the order of the approximation. References [17j, Il9|, 







20| include discussions of root search algorithms for the derivation of initial positions that 
reach the real axis at a given time. We will not describe all the details here but will just 
state the central idea. The complex root search exploits the assumption that the mapping 
x(0) 1-^ x{tf) is analytic. This property allows for an iterative process that detects the initial 
positions that correspond to real final positions. As demonstrated in references Q and 
in section UlI At for an arbitrary potential and final time, the mapping is only locally analytic. 
Generally, more than one initial position ends at a final position (whether real or complex). 
This makes the search for trajectories that end on the real axis more complicated but it has 
an important advantage in terms of interference effects. 

Our main observation is that the contribution of multiple trajectories in CTM can accu- 
mulate to an interference pattern. For simplicity we make the following assumption. Suppose 
that L trajectories end at final time t/ on real position x{tf) Then the final wavefunction is 
approximated by a superposition of contributions 

i^[x{tf),tf]^f2expi^^S'[xitf),tf]Y (2.29) 



where each trajectory (denoted by the index I) is associated with a phase S [x{tf),tf] 

N 

SHtf),tf] = J2fi'S^^[x{tf),tf], (2.30) 

that is calculated by the CTM equations of motion. In section IIIII we show that this as- 
sumption is too simplified and does not hold at all times and all positions. For example, 
for positions associated with a tunneling part of the wavefunction, only one of the multiple 
trajectories should be taken into account. A partial discussion on the superposition of con- 
tributions from complex trajectories appears in reference |l9j] in the GGWPD context. In a 
forthcoming paper [2^ we will explain an alternative derivation of the CTM in which the 
need to include multiple trajectories for certain times and positions becomes apparent 



III. NUMERICAL RESULTS 



In this section we examine numerically the CTM formulation allowing for the superposition 
of complex trajectories. For ready comparison the physical system we choose is identical 
to the one studied by Boiron and Lombardi (reference 17|] section IVB). The potential 
considered is a quartic double-well 



V{x) = 1.25 X W-\x^ 



The initial wavefunction is a Gaussian wavepacket 



400x' 



ip{x, 0) = exp 



-ao{x - XcY + -zPcix - Xc) + -7o 



(3.1) 



(3.2) 



where cto = 1, Xc = 0, pc = 5, 70 = — ^ In(^) and we take m = h = 1 (all quantities are 
given in atomic units). The initial conditions for the terms in the h power-expansion of the 
phase are 





;x,o) 


= iaoh{x — XcY + Pc{x — Xc) + 7o = ix'^ + 5x + 70, 


(3.3) 


d.So 


;x,o) 


= 2iaoh{x — Xc) + Pc = 2«x + 5, 


(3.4) 


dxxSo 


;x,o) 


= 2iaoh = 2i, 


(3.5) 


diSo 


;x,o) 


= 0, J > 3, 


(3.6) 




;x,o) 


= 0, J > 0, > 1, 


(3.7) 



where diSh = 



dxi ■ 



In section UlIAl we analyze the first order approximation of CTM (A^ = 1, 5* = 5*0 + ^5*1) 
and the properties of the trajectories. Section IIII Bl is dedicated to the next order of the 
approximation (A^ = 2, S = Sq + hSi + h?S2)- We omit an analysis of = since it is well 



presented in reference 



17| and only yields poor results. 



A. First Order approximation, A = 1 

The first order approximation of CTM requires the solution of eqs. fl2.15p . (12.161) . (12.171) . 
f l2.19p and fl2.2ip . The first two equations define the complex classical trajectories and the 
next three equations yield 5*0 and Si. We start by analyzing the complex classical trajec- 
tories. As mentioned above, the mapping a;(0) x{tf) is not one-to-one. For the quartic 
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potential, we found that three initial positions are mapped to every real final position at 
tf > 0. For short time scales this observation can be supported analytically. For general po- 
tentials or for longer time scales than we present here, more than three initial positions might 



19|, 



211 ] ■ In figures [T](a) and [T](b) we plot complex classical 



lead to the same final position 
trajectories for tf = 3 and tf = Q respectively. The initial positions of the trajectories can 
be divided into three groups referred to as 6ranc/ies[l^. One group of the initial positions 
is called the real branch and the other groups are called the secondary branches. The real 
branch is characterized by the property that it includes the initial position of a trajectory 
that propagates solely on the real axis. We refer to this trajectory as the real trajectory. 
It can be readily verified that for a Gaussian initial wavefunction there is only a single real 
trajectory that initiates at x(0) = Xc (see eqs.f l2.16]) . eq. fl2.15p and (13.41) ). In figJU^b) we de- 
pict the real trajectory explicitly. The secondary branches are defined simply as the groups 
of initial positions that do not belong to the real branch. Generally, the branches might be 
infinitely long curves in the complex plane. We will use the term branches to refer to the 
locus of initial positions that leads to final positions where the wavefunction is significantly 
different from zero. Hence, the branches are curves of finite length in the complex plane, 
although clearly there is some arbitrariness to their length. 

In fig{T](a) we see that at short time scales the secondary branches are centered far from 
neighborhood of the real axis. We can show analytically that for small times tf the initial 
positions that comprise the real branch obey |a;(0)| = 0{tf) whereas the secondary branches 
obey |a;(0)| = 0{jj). Note that the linear dependence of the initial momentum on position 
(eg. (13.41) ) allows trajectories with initial positions far from the real axis to reach a real final 
position in a short time. Unlike the secondary branches, the real branch is centered in the 
vicinity of the real axis at all times. The initial position x(0) = Xc is a fixed point of the 
real branch and prevents the real branch (recall that this is the locus of initial positions) 
from "straying" from the neighborhood of the real axis as the final time tf is increased. At 
intermediate times (time scales comparable to the time of the collision of the wavefunction 
with the barrier, 4 ^ tf ^ 7) secondary branch (1) reaches the vicinity of the real axis 
(figilt^b)) and at longer time scales it continues in the direction of the positive imaginary 
axis. As we demonstrate below, the proximity of secondary branch (1) to the real axis is 
closely related to the size of its contribution to the final form of the wavefunction and its 
role in interference effects. Secondary branch (2) does not reach the vicinity of the real axis 
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for any of the time scales specified below. The contribution of this branch to the absolute 
value of the final wavefunction (eg. (12.291) ) is negligible (in the order of 10~^^). Hence, from 
here on we ignore secondary branch (2) and refer to secondary branch (1) as the secondary 
branch. 

As we mentioned in section III D\ the existence of more than one branch motivates the 
attempt to superpose the contributions of the real branch and secondary branch in the final 
wavefunction 



tp[x{tf),tf]=i)u + i^s; t/'R = exp 



n 



exp 



n 



(3.8) 



where S'Reai and V'r are the phase and wavefunction associated with the real branch, and 
Ssec and ips correspond to the secondary branch. In figures [2t^a), [2t^b) and[2t^c) we compare 
the exact wavefunction with the numerical results obtained by applying CTM using a two- 
branch superposition. The figures indicate that when the wavefunction does not exhibit 
oscillations, the contribution of the real branch is sufficient to obtain a good approximation 
to the wavefunction. But at intermediate times, when the wavefunction exhibits interference 
effects, the contribution of both branches must be included. This last last observation applies 
in the spatial range up to the classical turning point (x ~ 24), beyond which the combined 
contribution diverges from the exact result. 

We turn to a closer inspection of this divergence. In figJHl we plot the contribution at 
tf = 6 of each individual branch and their superposition. Starting from the vicinity of 
X ^ 22, we observe an exponential increase of tps- For x ^ 23 we have a discontinuity of the 
approximation, as we discard the contribution of the secondary branch and include just the 



19j in the context of the 



real branch. A description of this divergence appears in reference 
GGWPD formulation. 

It is interesting to compare the time-dependence of the real and secondary branch con- 
tributions to the final approximation. A qualitative measure of the contribution of each 
branch is given by the imaginary part of the phase since 

S('S'Real) 



exp ( -^Real 



exp 



h 



(3.9) 



and a similar relation applies for ips and Q{Ssec)- In figureslD^a) andlD^b) we plot 9'(S'Reai) and 
^(5'sec) respectively for a series of final propagation times. We see that the secondary branch 
has a significant magnitude only at intermediate times. This observation coincides well with 
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FIG. 1: Complex classical trajectories with initial positions marked as circles and real final positions 
(marked as pluses) at (a) = 3 and (b) = 6. The trajectories arise from an initial Gaussian 
wavepacket propagating in a quartic double-well potential. The Gaussian is centered at x = 
and has positive initial momentum (the physical parameters are given in the text). In plot (a) we 
demonstrate that each final position arises from three initial positions. The initial positions are 
divided into a real branch and two secondary branches. The real branch is defined as incorporating 
a trajectory that remains on the real axis at all times. The real trajectory is specifically indicated 
in plot (b). 
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FIG. 2: A comparison between the exact quantum wavefunction and CTM (A^ = 1) with a two- 
branch superposition. The comparison is at a series of final propagation times specified by the 
numbers in the parentheses. The plots arise from an initial Gaussian wavepacket centered at x = 
with a positive average momentum, propagating in a quartic double-well potential (the parameters 
are given in the text), (a) Initially right-propagating wavefunction; (b) the reflected wavefunction; 
(c) a zoom on a section of (b). For = 5 in (a) and tf = 6 in (b) and (c) we plot the results both 
for just the real branch iV'Rp and for the combination of branches \ip^, + V'sp- The interference 
pattern obtained by superposing the contributions is clearly observed. 

the need to include the contribution of the secondary branch to the final wavefunction only 
at these times. The exponential growth of ips that is observed in figI3] is also apparent in 
figJHb), in the negative parts of the graphs for = 5 and tf = 6. The divergent magnitude 
of ips is in contrast to the finite magnitude of ip^, that is observed in fig|l](a). A discontinuity 
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FIG. 3: CTM approximation at = 6 for = 1. The contributions of each branch to the 
wavefunction are depicted by plotting [V'rI, jV'sl and |V'r + ■0s|- Note the exponential increase of 
ips begins around x ~ 22. For x ~ 23; we discard the contribution of the secondary branch and 
include just the real branch, leading to a discontinuity of the CTM approximation. 

in the derivative of 3(S'r) at tf = 5 and = 6 is also observed. This discontinuity appears 
slightly prior to the points where the contribution of the secondary branch begins to diverge. 

A close inspection of the complex trajectories at t/ = 5 and tf = 6, reveals an interesting 
property of the real trajectory: the real trajectory acts as a boundary between two "regimes" 
of complex trajectories comprising from the real branch. This can be seen in figJSl where the 
trajectories that initiate from Q'[x(0)] > are seen to reach the real x-axis at values lower 
than the real trajectory while trajectories with 53[x(0)] < seem to go past the barrier and 
reach the real x-axis at values higher than the real trajectory. These two regimes correspond 
to the two legs of the "v"-shaped graph of Q'(S'Reai) in figJHa): the trajectories arising from 
initial positions with S5[a;(0)] > correspond to the left leg of the "v" while trajectories with 
Q'[a;(0)] < correspond to the right leg of the "v". 



15 



B. Second Order approximation, N = 2 

In this section we analyze the effect of incorporating 5*2 in the CTM approximation. In 
addition to the five equations that are needed for obtaining the complex trajectories, 5*0 and 
^i, we need to solve eqs. (12.221) . (12.231) and (12.241) . In figE](a) we depict the approximate 
wavefunction for = 2 at t/ = 6. Comparing the N = 2 result with the = 1 result 
plotted in figlSl we conclude that other than an interval in the neighborhood of a; ~ 22.5, 
the N = 2 result (dashed line) lies on top of the exact result (solid line) and is significantly 
better then the = 1 result. For x ^ 23, where we incorporate solely the real branch 
contribution, the improvement in the approximation is graphically evident from the plots. 
For X < 22 we calculated the relative error between the absolute value of the approximations 
and the exact wavefunction using all the data points depicted in figsl3]and[6](a). The results 
are presented in figl6](b). For = 1 the mean relative error is 0.34% while for N = 2 the 
mean relative error is 0.11%. We see that the approximation worsens in the vicinity of the 
discontinuity of ipn. In the vicinity of x ~ 22.5, the N = 2 results are worse than the A^ = 1 
results; moreover, in the N = 2 case ips as well as t/'r exhibits a discontinuity. 

IV. SUMMARY 

In this paper we have presented a formulation of complex time-dependent WKB (CTDWKB) 
that allows the incorporation of interfering contributions to the wavefunction. The central 
idea in CTDWKB presented by Boiron and Lombardi[17] is to include both the amplitude 
and the phase in the lowest order term of the conventional time-dependent WKB method. 
The rationale behind this substitution is to treat the phase and the amplitude on equal 
footings in the limit h ^ 0. The benefits of the method are twofold. Firstly, CTDWKB 
exhibits accuracy superior to the conventional TDWKBp/7|]. Secondly, no singularities ap- 
pear in the integration of the equations of motion. The method has two main drawbacks. 
First, the trajectories that emerge obey the classical equations of motion but propagate in 
the complex plane (due to complex initial conditions), requiring analytic continuation of the 
quantum wavefunction. The second drawback is that the reconstruction of the wavefunction 
on the real axis requires a root search process. This process can be eased by exploiting the 
analytic mapping between initial and final position. 
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We have incorporated into the CTDWKB method the possibihty of contributions from 
multiple crossing trajectories. Boiron and Lombardi claim (section V in reference 171]) that 
they use the root search procedure "excluding de facto such double contributions" , although 
they appreciate the benefit that double contributions have in the GGWPD formulation. As 
we have demonstrated here, considering double contributions allows description of interfer- 
ence effects that are missing in the Boiron-Lombardi formulation of CTDWKB. Moreover, 
we have showed how to derive higher orders terms of the approximation in a straightforward 
manner. This process was applied for the derivation of a second order term in the CTDWKB 
approximation. The results for N = 2 were better than for = 1 except for a small interval 
in the vicinity of the classical turning point. It was also observed that even though there 
are no singularities in the integration of the CTDWKB equations of motion, a singularity 
appears in the real branch ipn at intermediate times. For = 2 an irregularity also appears 
in the part of the wavefunction associated with the secondary branch ips- We demonstrated 
that when a singularity appears in ipn (at intermediate times), the real trajectory acts as 
the boundary between two groups of trajectories associated with the real branch. Each of 
these groups contributes to a different side of the singularity. 

The CTDWKB formulation has several issues that require more comprehensive study. 
The most critical issue is to give an analytic explanation of the need to include the con- 
tributions from multiple classical trajectories (with zero relative phase) and why in some 
cases these contributions diverge. This will be dealt with in our forthcoming publication 



22| |. Some i nsig ht into the analytic structure of the complex classical trajectories was given 



in reference |l9j in the context of GGWPD; however, we believe that a more general under- 
standing of this structure is yet to be developed. This structure presumably is relevant to 
the question of when the CTDWKB formulation converges to the exact quantum mechani- 
cal result. We saw that in most parts of configuration space N = 2 performed better than 
A^ = 1, but in other parts of configuration space, where there were singularities, N = 2 
performed worse. What determines the position and time-dependence of these singularities 
in ipn at intermediate times? What is the relation between the singularities in CTDWKB vs. 
conventional time- dependent WKB? Is there any fundamental limitation on the time scale 
for which the method is accurate? Since WKB plays such a central role in quantum mechan- 
ics in general and in semiclassical mechanics in particular, we believe that these questions 
are of great general interest. The developments described in this paper together with the 
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answers to some of the above questions could make the time-dependent WKB formulation 
a competitive alternative to current time-dependent semiclassical methods. 

We wish to acknowledge David Kessler and Uzi Smilansky for useful discussions. This 
work was supported by the Israel Science Foundation (576/04). 
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FIG. 4: (a) 3'(S'Reai) and (b) 9(5scc) are depicted at a series of final propagation times (given in 
the parentheses). The results are limited to the spatial interval for which the absolute value of the 
exact wavefunction is significantly larger than zero. The imaginary part of the phase allows for a 
qualitative estimate of the contribution of each branch to the probability |^r + ^sPi see eq. (|3.9p . 
Figure (b) shows that $5(S'sec) drops below ~ 2 only for a finite interval of intermediate times. 
Therefore only for this range of times does the secondary branch makes a significant contribution 
to the wavefunction. 
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Real Branch, t =5 
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FIG. 5: Complex classical trajectories that correspond to the real branch at tj = 5. The real 
trajectory acts as a boundary between two "regimes" of the complex trajectories. Initial positions 
with Si[j;(0)] < seem to be go past the potential wall. The two "regimes" can be related to the 
singular behavior of the derivative of 9(5'ijea/) at intermediate times (figllja)). 
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t =6, N=2 




FIG. 6: (a) The second order (N = 2) CTM approximation is depicted for tj = 6. A discontinuity 
appears at x ~ 22.5 for both ipn and ips- (b) The relative error between the absolute value of the 
exact quantum wavefunction and the CTM approximation for N = 1 and N = 2, based on the 
data in figl3]and figEl^a). A comparison of the relative errors indicates a clear improvement when 
taking an additional order in the CTM approximation. 
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